

#include "glodef.h"
#include "errlog.h"
#include "f3io.h"
#include "fe.h"
#include "grid3d.h"

#include "glovar.h"

#ifdef USE_QGRID
#include "g3d_mesq.h"
#endif

#include "soah_common.h"

#include <stdlib.h>
#include <string.h>
#include <stdio.h>

#define SOAH_OUT_ROOT_DIR   "data"

static char _root_dir[FILENAME_MAX]=SOAH_OUT_ROOT_DIR;

const char *root_dir(void)
{
    return _root_dir;
}

static char _nodes_dir[FILENAME_MAX]=SOAH_OUT_ROOT_DIR "/nodes";

const char *nodes_dir(void)
{
    return _nodes_dir;
}

static char _logs_dir[FILENAME_MAX]=SOAH_OUT_ROOT_DIR "/logs";

const char *logs_dir(void)
{
    return _logs_dir;
}

static char _state_dir[FILENAME_MAX]=SOAH_OUT_ROOT_DIR "/state";

const char *state_dir(void)
{
    return _state_dir;
}

static char _adjoint_dir[FILENAME_MAX]=SOAH_OUT_ROOT_DIR "/adjoint";

const char *adjoint_dir(void)
{
    return _adjoint_dir;
}

void set_root_dir(const char *dir)
{
    size_t len=strlen(dir);
    if(len>FILENAME_MAX-10)
        error_msg("%s: path too long (%d) \n", __func__, (int)len);
    memcpy(_root_dir, dir, len+1);
    memcpy(_nodes_dir, dir, len+1);
    strcat(_nodes_dir, "/nodes");
    memcpy(_logs_dir, dir, len+1);
    strcat(_logs_dir, "/logs");
    memcpy(_state_dir, dir, len+1);
    strcat(_state_dir, "/state");
    memcpy(_adjoint_dir, dir, len+1);
    strcat(_adjoint_dir, "/adjoint");
}


void *solver_data;

const unsigned elemU=fe_bub;
const unsigned elemP=fe_p1;

double alpha=1.0;
double beta=1.0;

double total_cost;

/****************************************************************************/

void move_grid(void)
{

    if(flags.moving_grid==0)
        return;

    error_msg("moving_grid is disabled.\n");
    /*
        extern int dump_bin;
        if(dump_bin>0)
        {
            // copy current nodes' positions into old_nodes
            memcpy(alegrid->old_nodes, grid->nodes, grid->NN*sizeof(node_t));

            if(!load_nodes_time(nodes_dir()))
                error_msg("failed to load nodes\n");
            UpdateMeasures(grid);
        }
        else
        {
            move_alegrid_laplace(alegrid, TIME);

            grid->MAXH=GridSize(grid, grid->nodes);
            grid->MINQ=GridQuality(grid, grid->nodes);
    #if defined(QGRID_MEMBRANE)
            if((grid->MINQ>=1e-6)&&(grid->MINQ<QGRID_SMOOTH_THRESHOLD))
            {
                fork_QGrid(grid, grid->nodes, alegrid->can_move);
                grid->MAXH=GridSize(grid, grid->nodes);
                grid->MINQ=GridQuality(grid, grid->nodes);
                // we need to UpdateMeasures(grid) every time grid->nodes change
                UpdateMeasures(grid);
            }
    #endif
            ;
            if(grid->MINQ<1e-6)
                error_msg("Grid quality %g is too low.\n", grid->MINQ);
            PrintLog("(INFO) current grid: size %g     quality %g\n",
                     (double) grid->MAXH, (double) grid->MINQ);

        }

        // compute the norm of the gradient of the mesh velocity
        vector nd=get_ndisp(alegrid);
        scale_vec(nd, 1.0/TStep);

        node_t *nds=grid->nodes;
        grid->nodes=alegrid->ref_nodes;
        UpdateMeasures(grid);
        double norm=norm_grad_vec(grid, fe_p1, nd);
        PrintLogNorm("grad MeshVel", "L2", norm);
        grid->nodes=nds;
        UpdateMeasures(grid);

        if(massP) fe_matr_assemble(grid, 1, &massP);
        if(massU) fe_matr_assemble(grid, 1, &massU);
     */
}

/****************************************************************************/

double CostFunctional(void)
/* compute value of cost functional at this time */
{
    double cost;
    double NGU, NCU, NDU;
    NGU=norm_grad_vec(grid, elemU, u);
    NDU=norm_div_vec(grid, elemU, u);
    NCU=norm_curl_vec(grid, elemU, u);
    PrintLogNorm("div  U", "", NDU);
    PrintLogNorm("curl U", "", NCU);
    PrintLogNorm("grad U", "", NGU);
#if (COST_FUNC_SQUARE == 0)
    cost=sqrt(alpha*sqr(NGU)+beta*sqr(NCU));
#else
    cost=0.5*(alpha*sqr(NGU)+beta*sqr(NCU));
#endif
    PrintLog("       %10s   %-20s = %e\n", "cost", "", cost);
    return cost;
}

/****************************************************************************/

void set_defaults(void)
{
    /* set default flags */
    flags.advection=adv_semi;
    flags.adj_advection=adv_semi;
    flags.time_order=1;
    flags.transfer=transfer_none;
    flags.moving_grid=FALSE;
    flags.stationary=FALSE;
    flags.convergence=FALSE;
    flags.conv_adjoint=FALSE;
    flags.periodic=FALSE;
    flags.shape_opt=FALSE;
    flags.need_adjoint=FALSE;
    flags.solver=solver_cs;
    flags.use_ale=TRUE;

    /* set default parameter values */
    Re=1.0;
    Fr=0.0; /* No gravity */
}

unsigned on_ival=0,
        on_oval=0,
        on_memb=0,
        on_pipe=0,
        on_chbr=0;

void soah_check_heart_geometry(void)
{
    if(on_ival==0||on_oval==0||on_memb==0||
        (on_ival&on_oval)!=0||(on_ival&on_memb)!=0||
        (on_ival&on_pipe)!=0||(on_ival&on_chbr)!=0||
        (on_oval&on_memb)!=0||(on_oval&on_pipe)!=0||
        (on_oval&on_chbr)!=0||(on_memb&on_pipe)!=0||
        (on_memb&on_chbr)!=0||(on_pipe&on_chbr)!=0)
        error_msg("%s: must specify valid heart-geometry.\n", __func__);
    if(on_chbr==0)
        on_chbr= ~(on_ival|on_oval|on_memb|on_pipe);
    else
        if((on_chbr&on_ival)!=0||(on_chbr&on_oval)!=0||
        (on_chbr&on_memb)!=0||(on_chbr&on_pipe)!=0)
        error_msg("%s: chamber overlaps with other boundaries.\n", __func__);
}
